Species change at pixel level

library(tidyverse)
library(terra)
library(sf)
library(rnaturalearth)
library(here)
library(iucnredlist) ### remotes::install_github("IUCN-UK/iucnredlist")
library(janitor)
library(cowplot)

See how species change at the pixel level

Will create a map of Soresen’s Similiarity Index at the pixel level

# Functions 
source(here("common_fxns.R"))

# Read in Lawler data
spp_mapfiles <- list.files(here_anx('/data_lawler_2020/species_projections'), 
                           full.names = TRUE, recursive = TRUE, pattern = '.tif')
spp_all_df <- data.frame(f = spp_mapfiles,
                         spp = basename(spp_mapfiles) %>%
                           str_remove('(_gf85|_in85|_mc85|_bias).+') %>%
                           str_replace_all('_', ' '),
                         tax = basename(dirname(spp_mapfiles)))

# Read in AVONET dataset (reading in BirdLife tab, which should match with IUCN names)
avonet <- read_csv(here("data-laa", "AVONET-birdlife-tab-supplementary-dataset-1.csv")) %>% 
  clean_names()

# Filter for just Lawler birds
birds_all_df <- spp_all_df %>% 
  filter(tax == "Birds")

birds_all <- birds_all_df %>% 
  group_by(spp, tax) %>% 
  summarize(files = n())

Gather maps and aggregate

Following the method in 2_combine_spp_scenarios.qmd: For all the Lawler birds species (301) gather all future scenario maps and aggregate, summing and dividing by 3 to make a “probability” map for future distributions, and crop to just the lower 48 US states (and a little Canada and Mexico). Then crop the historical maps to the same extent.

Aggregate future projections

For each species in the sample, aggregate the three future scenarios into a single map. Crop to the lower continental US states.

### set up extent for continental US
us_extent <- ext(c(xmin = -130, xmax = -65, ymin = 23, ymax = 50))

### Aggregate over the vector of species names
spp_vec <- birds_all_df$spp %>% unique()

future_rast <- lapply(spp_vec, FUN = agg_scenario, df = birds_all_df) %>%
  setNames(spp_vec) %>%
  rast()

x <- future_rast[[1]]
#plot(x, main = names(future_rast)[1])

crop historical maps

hist_rast <- lapply(spp_vec, FUN = agg_scenario, df = birds_all_df, 
                    scenario = 'historical') %>%
  setNames(spp_vec) %>%
  rast()

x <- hist_rast[[1]]
#plot(x, main = names(hist_rast)[1])

Generate jurisdiction map

Let’s generate a simple jurisdiction map - US states, and Canada and Mexico.

states_sf <- ne_states(returnclass = 'sf', country = c('Mexico', 'Canada', 'United States of America')) 

state_id_df <- states_sf %>%
  st_drop_geometry() %>%
  mutate(juris = ifelse(adm0_a3 == 'USA', name, admin)) %>%
  dplyr::select(id = gn_id, juris)
### plot(states_sf %>% select(iso_a2))
states_rast <- rasterize(states_sf, hist_rast, field = 'gn_id')
#plot(states_rast)

# Function to get pixels instead of jurisdictions
extract_jurisdictions_pixels <- function(spp, future_r, hist_r, juris_r) {
  ### spp <- spp_vec[1]
  ### future_r = future_rast; hist_r = hist_rast; juris_r = states_rast
  s_f_r <- future_r[[names(future_r) == spp]]
  s_h_r <- hist_r[[names(hist_r) == spp]]
  s_df <- data.frame(values(s_h_r),
                     values(s_f_r),
                     values(juris_r)) %>%
    setNames(c('historical', 'future', 'id')) %>%
    mutate(pixel_id = row_number()) %>% 
    filter(!is.na(id)) %>%
    filter(!is.na(historical) | !is.na(future)) #%>%
    #group_by(id) %>%
    #summarize(n_historical = sum(historical > .25),
              #n_future = sum(future > .5))

  return(s_df)
}

Rough pass at classification?

For each spp let’s create a dataframe of jurisdictions and cell values…

state_spp_list <- lapply(spp_vec, FUN = extract_jurisdictions_pixels,
                         future_r = future_rast, 
                         hist_r = hist_rast,
                         juris_r = states_rast) %>%
  setNames(spp_vec) %>%
  bind_rows(.id = 'spp') %>%
  left_join(state_id_df, by = 'id') #%>%
  ### note, Canada and Mexico have multiple jurisdictions (provinces and states)
  ### so sum up all those instances
  #group_by(spp, juris) %>%
  #summarize(n_historical = sum(n_historical),
            #n_future = sum(n_future), 
            #.groups = 'drop')

Now that we have each species at the pixel level (pixel id), how many species are in each pixel for historical and future?

pixel_spp <- state_spp_list %>% 
  group_by(pixel_id, juris) %>% 
  summarize(n_spp_historical = sum(historical > .25),
            n_spp_future = sum(future > .5)) %>% 
  mutate(n_spp_diff = n_spp_future - n_spp_historical)

# Vast majority of pixels are losing total # of species (simply from numbers perspective)
ggplot(data = pixel_spp) +
  geom_histogram(aes(x = n_spp_diff),
                 bins = 50) +
  labs(title = "Bird species gain from historic to future by pixel",
       x = "species change")

Next step - sorensen similarity index (test a few pixels)

# Test out with a single pixel (141325)
pxl_141325 <- state_spp_list %>% 
  filter(pixel_id == 141325) %>% 
  # Assign future pixels > 0.5 to spp prescence 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0))

pxl_si <- pxl_141325 %>% 
  group_by(pixel_id) %>% 
  summarize(shared = sum(historical == 1 & future == 1),
            historical_spp = sum(historical),
            future_spp = sum(future),
            ssi = (2 *shared)/(historical_spp + future_spp))

# Try all pixels
pxl_ssi <- state_spp_list %>% 
  # Assign future pixels > 0.5 to spp prescence 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0)) %>% 
  group_by(pixel_id, juris) %>% 
  summarize(shared = sum(historical == 1 & future == 1),
            historical_spp = sum(historical),
            future_spp = sum(future),
            ssi = (2 *shared)/(historical_spp + future_spp))

# Distribution of Soresen Similiarity Index for bird species
# Mean is 0.40 (moderate similarity between historic and future)
ggplot(data = pxl_ssi) +
  geom_histogram(aes(x = ssi),
                 bins = 50) +
  geom_vline(xintercept = mean(pxl_ssi$ssi)) +
  annotate("text", label = "mean = 0.40", x = 0.30, y = 15000) +
  #facet_wrap(~juris) +
  labs(title = "Soresen Similarity Index btwn historic & future by pixel (n = 152,594)",
       x = "SSI")

# SSI stats by jurisdiction
ssi_juris <- pxl_ssi %>% 
  group_by(juris) %>% 
  summarize(mean_ssi = mean(ssi),
            median_ssi = median(ssi),
            max_ssi = max(ssi),
            min_ssi = min(ssi)) %>% 
  ungroup()

ssi_juris %>% 
  mutate(juris = fct_reorder(juris, mean_ssi, .desc = TRUE)) %>%
  ggplot() +
  geom_col(aes(x = juris, y = mean_ssi)) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  labs(y = "mean pixel SSI",
       x = "jurisdiction",
       title = "Bird species mean pixel SSI by jurisdiction")

Create sf of raster ids to map the SSI by pixel

# Create a blank raster with same extent as x = hist_rast[[1]] 
id_rast <- rast(x)

# Assign IDs and then convert to sf
values(id_rast ) <- 1:ncell(x)

id_rast_sf <- id_rast %>% 
  stars::st_as_stars() %>% 
  st_as_sf() %>% 
  rename("pixel_id" = "Acrocephalus familiaris")

# Merge in the SSI
id_rast_sf_ssi <- inner_join(id_rast_sf, pxl_ssi)

ggplot() +
  geom_sf(data = id_rast_sf_ssi %>% 
            filter(juris == "California"), aes(geometry = geometry, fill = ssi, color = ssi), size = 0) + 
  scale_fill_viridis_c(limits = c(0,1)) + 
  scale_color_viridis_c(limits = c(0,1)) +
  labs(title = "Birds in CA Sorenen's Similiarity Index")

ggplot() +
  geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, color = ssi, fill = ssi)) +
  scale_color_viridis_c() + 
  scale_fill_viridis_c() + 
  labs(title = "Birds Sorenen's Similiarity Index")

In general, SSI values are as follows (https://www.statology.org/how-to-calculate-and-interpret-the-sorensen-dice-index/):

• 0.80-1.00: Very high similarity (communities share most species)
• 0.60-0.79: High similarity (substantial species overlap)
• 0.40-0.59: Moderate similarity (some shared species)
• 0.20-0.39: Low similarity (few shared species)
• 0.00-0.19: Very low similarity (minimal species overlap)

historic_spp_num <- ggplot() +
  geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, fill = historical_spp, color = historical_spp)) + 
  scale_fill_viridis_c(option = "plasma", limits = c(1,180)) + 
  scale_color_viridis_c(option = "plasma", limits = c(1,180)) +
  geom_sf(data = id_rast_sf_ssi %>% filter(historical_spp == 0), aes(geometry = geometry), fill = "gray40", color = "grey40") + 
  labs(title = "Historic number of bird species",
       subtitle = "(0 species in gray)",
       fill = "count",
       color = "count")

future_spp_num <- ggplot() +
  geom_sf(data = id_rast_sf_ssi, aes(geometry = geometry, fill = future_spp, color = future_spp)) + 
  scale_fill_viridis_c(option = "plasma", limits = c(1,180)) + 
  scale_color_viridis_c(option = "plasma", limits = c(1,180)) +
  geom_sf(data = id_rast_sf_ssi %>% filter(future_spp == 0), aes(geometry = geometry), fill = "gray40", color = "grey40") + 
  labs(title = "Future number of bird species",
       subtitle = "(0 species in gray)",
       fill = "count",
       color = "count")

plot_grid(historic_spp_num + theme(legend.position = 'bottom'), 
          future_spp_num + theme(legend.position = 'bottom'))

ggplot() +
  geom_sf(data = id_rast_sf_ssi %>% 
            filter(juris == "Canada"), aes(geometry = geometry, fill = ssi, color = ssi), size = 0) + 
  scale_fill_viridis_c(limits = c(0,1)) + 
  scale_color_viridis_c(limits = c(0,1)) +
  labs(title = "Birds in Canada Sorenen's Similiarity Index")

Add in Avonet data and investigate trait change at pixel level

Let’s start with the birds that easily merge with the Avonet dataset (270)

# See how well the Avonet trait data merges in 
# 270 of 301 have matches - could work with these species for now since they should be able to match with IUCN
birds_avonet <- left_join(birds_all, (avonet %>% rename(spp = species1)))

lawler_avvonet_birds <- birds_avonet %>% 
  drop_na(sequence)

lawler_avvonet_birds_vec <- lawler_avvonet_birds$spp

Now, will add in the bird species that did not merge as easily for a total of 299

# Read in crosswalk of lawler to avonet species (created with input from gemini)
lawler_avonet_xwlk <- read_csv(here("data-laa", "lawler_avonet_xwlk.csv")) %>% 
  clean_names() %>% 
  dplyr::select(lawler_spp:name_status) %>% 
  drop_na()

# Match the lawler avonet xlk with trait data
# The two extinct species drop out
lawler_avonet_traits_missing <- inner_join(lawler_avonet_xwlk, 
                                   avonet %>% rename(avonet_spp = species1))

# Bind the easily matched birds and the birds that were missing in the merge
lawler_avonet_all <- rbind(
  lawler_avvonet_birds %>% dplyr::select(-tax, -files),
  lawler_avonet_traits_missing %>% dplyr::select(-avonet_spp, -name_status)
)

pxl_spp_list_trait <- inner_join(state_spp_list, lawler_avonet_all) %>% 
    filter(!is.na(id)) %>%
    filter(!is.na(historical) | !is.na(future))

Create map

us_spp_trait <- pxl_spp_list_trait %>% 
  dplyr::select(spp:juris, habitat:primary_lifestyle) %>% 
  filter(historical > 0.25 | future > 0.5) %>% 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0)) %>% 
  mutate(habitat_density = case_when(habitat_density == 1 ~ "Dense",
                                     habitat_density == 2 ~ "Semi-open",
                                     habitat_density == 3 ~ "Open")) %>% 
  mutate(migration = case_when(migration == 1 ~ "Sedentary",
                               migration == 2 ~ "Partial",
                               migration == 3 ~ "Migratory"))

# Combine the categorical columns into one group column to create unique groups
# Count up these trait groups by pixel
us_spp_trait_groups <- us_spp_trait %>% 
  unite(trait_group, migration, habitat_density, habitat, trophic_level, trophic_niche, primary_lifestyle) %>% 
  group_by(pixel_id, juris, trait_group) %>% 
  summarize(historical_trait_group_count = sum(historical),
            future_trait_group_count = sum(future)) %>% 
  ungroup()

# If we apply the sorensen index to traits, then we care about the presence / absence of a trait group and not the # of species
us_spp_trait_groups_ssi <- us_spp_trait_groups %>% 
  mutate(historical_trait_group_presence = case_when(historical_trait_group_count == 0 ~ 0,
                                                     historical_trait_group_count != 0 ~ 1),
         future_trait_group_presence = case_when(future_trait_group_count == 0 ~ 0,
                                                     future_trait_group_count != 0 ~ 1)) %>% 
  group_by(pixel_id, juris) %>% 
  summarize(shared = sum(historical_trait_group_presence == 1 & future_trait_group_presence == 1),
            historical_trait_groups= sum(historical_trait_group_presence),
            future_trait_groups = sum(future_trait_group_presence),
            trait_group_ssi = (2 *shared)/(historical_trait_groups + future_trait_groups)) %>% 
  ungroup()

# Merge in pixel geography for map
id_rast_sf_us_ssi <- inner_join(id_rast_sf, us_spp_trait_groups_ssi)

# Create map (299 bird species)
ggplot() +
  geom_sf(data = id_rast_sf_us_ssi, aes(geometry = geometry, fill = trait_group_ssi, color = trait_group_ssi), size = 0) + 
  scale_fill_viridis_c(limits = c(0,1)) + 
  scale_color_viridis_c(limits = c(0,1)) +
  labs(title = "Bird trait groups Sorenen's Similiarity Index")

What’s the change from species SSI to trait SSI

# Filter the species dataframe down to on the avonet birds
pxl_ssi_avonet_spp <- state_spp_list %>% 
  # Assign future pixels > 0.5 to spp prescence 
  filter(spp %in% lawler_avvonet_birds_vec) %>% 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0)) %>% 
  group_by(pixel_id, juris) %>% 
  summarize(shared = sum(historical == 1 & future == 1),
            historical_spp = sum(historical),
            future_spp = sum(future),
            ssi = (2 *shared)/(historical_spp + future_spp))

# Merge in the SSI
id_rast_sf_ssi_avonet_spp <- inner_join(id_rast_sf, pxl_ssi)

Would expect trait SSI to generally to be greater than species SSI

pxl_ssi_dif <- left_join(us_spp_trait_groups_ssi %>% 
                           dplyr::select(pixel_id, trait_group_ssi),
                         pxl_ssi %>% 
                           dplyr::select(pixel_id, ssi)) %>% 
  mutate(ssi_diff = trait_group_ssi - ssi)

# Merge in the pixel geography
id_rast_sf_us_ssi_diff <- inner_join(id_rast_sf, pxl_ssi_dif)

ggplot() +
  geom_sf(data = id_rast_sf_us_ssi_diff, aes(geometry = geometry, fill = ssi_diff, color = ssi_diff), size = 0) + 
  scale_color_distiller(palette = "BrBG", direction = 1, limits = c(-0.4, 0.4)) + 
  scale_fill_distiller(palette = "BrBG", direction = 1, limits = c(-0.4, 0.4)) +
  # Show areas 0 species and trait overlap in pink
  geom_sf(data = id_rast_sf_us_ssi_diff %>% 
            filter(ssi == 0 & trait_group_ssi == 0), 
          aes(geometry = geometry), fill = "pink", color = "pink") +
  # Show areas with full species and trait overlap in navy
  geom_sf(data = id_rast_sf_us_ssi_diff %>% 
            filter(ssi == 1 & trait_group_ssi == 1), 
          aes(geometry = geometry), fill = "navy", color = "navy") +
  labs(title = "Trait SSI - Species SSI",
       subtitle = "(pink is 0 trait and 0 species SSI; navy is 1 trait and 1 species SSI)")

#Let’slook closer at the pixels where the difference is 0 Zeros may be areas where there are 0 species SSI and 0 trait SSI (no similarity btwn historic and future), 1 species SSI and 1 trait SSI (complete similarity btwn historic and future), or nonzero values.

pxl_ssi_diff_0 <- pxl_ssi_dif %>% 
  filter(ssi_diff == 0)

pxl_ssi_diff_0_summary <- pxl_ssi_diff_0 %>% 
  mutate(diff_type = case_when(trait_group_ssi == 0 & ssi == 0 ~ "zero",
                               trait_group_ssi == 1 & ssi == 1 ~ "one",
                               TRUE ~ "non-zero")) %>% 
  group_by(diff_type) %>% 
  summarize(count = n()) %>% 
  ungroup()

zero_trait_species_ssi <- pxl_ssi_diff_0 %>% 
  filter(trait_group_ssi == 0 & ssi == 0 )

zero_trait_species_ssi_pxls <- zero_trait_species_ssi$pixel_id

# Look at the zero ssi pixels

zero_species_ssi <-state_spp_list %>% 
  # Assign future pixels > 0.5 to spp prescence 
 # filter(spp %in% lawler_avvonet_birds_vec) %>% 
  filter(pixel_id %in% zero_trait_species_ssi_pxls) %>% 
  mutate(future = case_when(future > 0.5 ~ 1,
                           TRUE ~ 0)) %>% 
  # Assigning presence / absence for these birds
  filter(historical == 1 | future == 1)

Plotting a few of the North Dakota birds

# Seems like it's one of the few birds showing up in future North Dakota (even though it seems unlikely to be there based on research)
plot(hist_rast$`Xenus cinereus`, main = "Xenus cinereus - historical")

plot(future_rast$`Xenus cinereus`, main = "Xenus cinereus - future")

plot(hist_rast$`Parus hudsonicus`, main = "Parus hudsonicus - historical")

plot(future_rast$`Parus hudsonicus`, main = "Parus hudsonicus - future")

plot(hist_rast$`Anser rossii`, main = "Anser rossii - historical")

plot(future_rast$`Anser rossii`, main = "Anser rossii - future")

plot(hist_rast$`Tympanuchus pallidicinctus`, main = "Tympanuchus pallidicinctus - historical")

plot(future_rast$`Tympanuchus pallidicinctus`, main = "Tympanuchus pallidicinctus - future")